Drying Patterns: Sensitivity to Residual Stresses 
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Volume alteration in solid materials is a common cause of material failure. Here we investigate 
the crack formation in thin elastic layers attached to a substrate. We show that small variations in 
the volume contraction and substrate restraint can produce widely different crack patterns ranging 
from spirals to complex hierarchical networks. The networks are formed when there is no prevailing 
gradient in material contraction whereas spirals are formed in the presence of a radial gradient in 
the contraction of a thin elastic layer. 



Introduction. Desiccation is known to produce com- 
plex networks of shrinkage-cracks in starch-water mix- 
tures or clays [l|, [3, 0, 0]- In concrete small cracks are 
often formed by the preparatory drying process and by 
the later ingress of reactive reagents. Similarly in na- 
ture, the infiltration of fluids and chemical reagents into 
rocks generate internal stresses that form intricate pat- 
terns of pervasive cracks [5]. Typically the stress is gen- 
erated from local volume changes. Fractures are also ob- 
served in thin films attached to a substrate. Experiment 
on films have revealed intricate patterns ranging from 
the hierarchical structure typically observed in mud and 
concrete to spiral shaped cracks Q, @. In spin-coating 
a fluid droplet is added at the center of a rotating sub- 
strate and is spread by centrifugal forces to cover the 
full substrate. During the drying and curing of the sys- 
tem, chemical bonds are formed between the coating and 




Iff** 



FIG. 1: Color online. Four different stages in the evolution of 
a hierarchical crack network. The total contraction was 9% 
and cracks were nucleated inside domains whenever the max- 
imum principal stress exceeded a c = 0.85 in arbitrary units. 
The substrate restraining force was drawn from a uniform 
distribution with 10% disorder and a mean of unity. 



the substrate. In this process the coating often shrinks 
and tensile stresses are produced [7[ and if one is less 
careful, the stress may result in an unwanted cracking 
of the coating. Due to the spinning of the system, the 
residual stress of the coating may also contain inherent 
shear components. Other mechanisms like anisotropic 
drying rates can also leave behind remnant shear. In 
cases where the thickness of the drying specimen is small 
and the contraction fairly uniform, i.e. no residual shear 
stresses, the growing cracks typically form an intricate hi- 
erarchical pattern. The pattern is the result of a cascade 
of successive cracks (which is supported by experimen- 
tal evidences [2]); at each fragmentation stage, a crack is 
forming which divides a mother fragment of area A into 
two daughter fragments of areas A\ and A2 respectively, 
with area conservation (A = A\ + A2). It is worth notic- 
ing that in principle the trajectory of the crack that is 
dividing the mother fragment can be anything, and is de- 
termined from the shape and size of the mother domain 
and the inherent material disorder. 

In order to analyze the spiral and hierarchical cracks, 
we consider a system consisting of a thin elastic layer 
attached to an elastic substrate. Under plane stress con- 
ditions we have that the in-plane strain tensor e^- in the 
layer is related to the stress tensor <j^- by 
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where E is Young's modulus and (3 (in the absence of 
external stress) is a measure of the free volume change 
caused by e.g. drying or thermal expansion of the thin 
elastic layer. Whenever the film is displaced from its 
equilibrium position by a local displacement u the elas- 
tic substrate tries to restore the film by a force f(u). 
For small displacements we assume that this force is lin- 
early proportional to — u. In general it is assumed that 
volume alteration in the film happens on a time scale 
much larger than the time required for elastic waves to 
propagate across the system and the system is therefore 
assumed to always be in elastostatic equilibrium. The 
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force balance therefore assumes the form 
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where \i is the constant of proportionality of the substrate 
restoring force. For small deformations the strain follows 
from the displacement via the relations = (djUi + 
diUj)/2. Combining this relation with the force balance 
Eq. ([2]) and stress-strain relations Eq. (pQ) we achieve the 
following equation for the displacement 
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We now provide an estimate of the typical stress en- 
countered during volume alteration of thin films with a 
linear spatial extend of size R. To that end, we shall con- 
sider the maximum stress for a circular domain of radius 
R located at the center of coordinates and with vanishing 
stress at the boundaries. The displacement field is for a 
uniform material contraction (3 found as a solution to the 
radial symmetric version of Eq. ([3|) 
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where the material specific constant a is given by a = 
(1 - v 2 )n/E. Multiplying both sides of Eq. (|4]) by r 2 and 
rescaling r with ^/a yields the modified bessel differential 
equation. The solution for o rr (R) = and u r (0) = is 
given by 
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Here J n , n = 0, 1 is the modified bessel function of the 
first kind. From the displacement field, the stress in 
cylindrical coordinates follows from the expressions 
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Note that by the symmetry of the problem the shear 
stress vanishes. The breaking of this symmetry will be 
important for the formation of the spiral crack patterns 
presented below. The stress components have their max- 
imum (absolute value) at the middle of the circular do- 
main and are given by 
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The magnitude of the stress components monotonically 
increases with R and in the limit R — » 00 the stress com- 
ponents achieve the value —E(3/(l — v). In the limit 
R = the stress becomes —E/3/2. From Eq. (0) we can 



now provide an estimate of a critical domain size that will 
fracture under a predefined yield stress. As long as this 
yield stress is lower than the material stress fracture will 
form and grow. Depending on the material contraction 
the fractures may develop into spiral shaped patterns or 
hierarchical networks. In both cases the maximum stress 
is reduced by the propagating crack and only when the 
stress drops below the yield stress the fracturing stops. In 
Fig. [I] we show a fracture network resulting from numer- 
ical solutions as explained below, formed from an initial 
contraction of the substrate and a predefined yield stress 
level. According to Eq. ([7j) each domain division reduces 
the stress and an average linear size (R) of the domains 
can be found for a given yield stress by inverting Eq.([7j). 

The non-uniformity of the elastic layer and the com- 
plex boundary conditions make it hard and often impos- 
sible to find an analytical solution to the displacement 
equation. Therefore we have implemented a numerical 
method based on the Galerkin finite element discretiza- 
tion using an adaptive triangular meshing. In the vicinity 
of a propagating crack tip we highly increase the resolu- 
tion by decreasing locally the area of the triangular ele- 
ments and thereby allow for an accurate computation of 
the stress intensity factors of the propagating crack. The 
drying process is simulated by applying a body force to 
the elements, i.e. we shift the equilibrium position by 
adding an extra force term on the right hand side of Eq. 
([2j) . In that way we can readily add disorder into the sys- 
tem by selecting the magnitude of the local body force 
from a random distribution. In the simulations on hi- 
erarchical fracture networks presented below, we use a 
uniform distribution with unit mean. 

Crack Initiation and propagation. Here we 
present in detail how we nucleate cracks and model their 
evolution. First we find the points which have the high- 
est stress and exceed the critical value. The stress is 
determined along the principal axes of the stress matrix 
<Jij, i.e. the principal stress. Whenever the yield stress 
is exceeded, we nucleate at the point of yielding a small 
semi elliptical void with an eccentricity of 0.998. The 
major axis of the ellipse is aligned in the direction of the 
maximum principal stress. We allow the crack to evolve 
according to the Griffith criterion and the principle of 
local symmetry, i.e. the crack will grow in a direction 
such as to annul the local shear component at the crack 
tip. At each step of propagation we compute the stress 
in every element near the crack tip and find the stress 
intensity factors from a best fit to the equations [3] • 
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Here r, 6 are local polar coordinates with respect to the 
crack tip with measured from the line following the 
direction of the crack. 099 and o r 9 are the circumferential 
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FIG. 2: Color online. The distribution density of \\ — 1/2 1 
for simulations with 5%, 10% and 15% disorder on the sub- 
strate restraining. The distributions are averaged over do- 
mains formed at the 6th generation of cracks. No significant 
variation is seen at these fairly low levels of disorder. The 
black line on top represents a best fit with an exponential 
distribution exp( — \x — l/2\/a) where a = 0.03. In the inset 
we show for the same data the cumulative distributions of the 
domain areas. The dashed line on top is an estimate of the 
distributions considering the individual domain divisions to 
be uncorrelated (see text). 



tensile stress and the shear stress, respectively. Kj and 
Kji are the unknown stress intensity factors for mode 
I and II, respectively. The principle of local symmetry 
is satisfied if the crack grows in a direction given by an 
angle a where Ku — > 0. Suppose that the crack forms 
an infinitesimal kink at an angle a from the old direction 
of the crack, we can define the local mode I and mode II 
stress intensity factors, 
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Whether the crack propagates or not is dictated by the 
Griffith criterion, i.e. the energy balance of the energy 
release rate into the crack tip region must balance the dis- 
sipation involved in the crack propagation. For a kinked 
crack, the energy release rate is, 

G(a) = (KUa)+Kh(a))/E . (11) 
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FIG. 3: Simulation of spiral cracks for various values of the 
material contraction /3 and the rotation angles. In all the 
panels, the crack was initiated at the center and was allowed 
to propagate until it reached the outer boundary. Note that 
the smaller the shear stress is the more pronounced is the 
spiralling. 



the maximum of the strain energy release rate 
dG(a)/da=0 is equivalent to Ku(a)=0 or dKi(a)/da = 
0, thus the new direction of the crack ao corresponds 
to the point where Kj(ao) exhibits a maximum and 
K H (a ) = [13]. Applying the latter to Eq. (10]) yields, 
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We emphasize that in this model both the cracking time 
and the area of the fragmented elements depend only 
on material contraction and initial disorder. In a natu- 
ral system this may not always be the case since mate- 
rial properties and disorder can evolve in time. Uniform 
contraction of the elastic layer produces homogenous hi- 
erarchical crack patterns. One can in this case argue 
that the effect of the crack is to partition the mother 
area A (of generic shape) into two areas A\ = \A and 
A2 = (1 — x)A, where < x < 1 is a random variable 
whose distribution (that must be symmetric under the 
transformation \ l— ^ 1 — x) is unknown. In Fig. [2] is the 
distribution of |% — 1/2 1 = \A\ — A^jA shown together 
with a best fit to an exponential distribution. Although 
the domain areas are correlated to their mother domains, 
the exponential distribution of \ allows for a simple es- 
timate of the area distribution by neglecting the correla- 
tion. That is the areas at the n'th generation level are can 
be determined by a product of n random numbers drawn 
from the exponential distribution, i.e. A^ = Ao Yl™ 
1 In Fig. [2j we show in the inset a distribution for domain 
areas at the 6th generation together with the distribution 
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FIG. 4: Color online. The figure shows the relation between k 
and the rotation angle for four different values of the material 
contraction /3. n is computed as the exponent of a best fit to 
a logarithmic spiral, r(0) = roe K< ^ . The rotation angle is the 
prefactor 60 used in the expression for the relative rotation 
between the substrate and the thin film. The inset shows a 
data collapse of k/3 for the same curves and is in agreement 
with the simple scaling form k ~ 9q//3. 



of the product of random numbers. After a few number of 
generations this distribution will approach a log-normal 
distribution. The deviation from the exponential distri- 
bution for larger values of \x — 1/2 1 will for an increasing 
number of fracture generations lead to a less good fit 
using the exponential distribution as an approximation. 

Spirals. We now investigate what happens when the 
contraction is non-uniform and has smooth gradients. 
Experiments on thin films attached to a substrate by a 
spin-coating technique reveal a broad range of crack pat- 
terns @ ranging from networks of cracks to single cracks 
spiraling outwards from their site of nucleation. The nu- 
cleation is usually taking place at localized sites with high 
stress typically generated from small defects or inclusions 
in the material. If the material contraction is uniform 
in the neighborhood around the crack, the crack would 
propagate straight towards the material boundary where 
it would curve to meet the boundary at a right angle. 
However if the material contraction increases smoothly 
away from the site of nucleation, straight cracks would 
be unstable and would start to curve. In that way circu- 
lar shaped cracks can be formed. During the preparation 
of thin film coatings (such as spin-coating) it is not un- 
common to have a minor residual shear stress. The shear 
stress breaks the symmetry of the system and the circular 
crack may then turn into a spiral. If we alter the con- 
traction such that it increases linearly away from a given 
site of crack nucleation, (e.g. follows a simple linear form 
(3(r) = /3r, and add a small shear stress by rotating the 



elastic layer relative to the underlying substrate with an 
angle eq (r) — #or 7 ), the crack would propagate along 
a spiral trajectory. Different powers of 7 result in the 
formation of different spirals. For the simulation of the 
crack propagation we use an initially circular symmetric 
system. A crack is then initiated at the center of the cir- 
cle and is allowed to propagate according to the Griffith 
criterion until it reaches the boundary. The results using 
a rotation of the substrate by a power 7=1/2 are shown 
in Fig. [3] using various values for the pref actors #0 and /?, 
respectively. The cracks have a shape that fit well a log- 
arithmic spiral, i.e. they have a form r(0) — Aexp(ft#) 
where k depends on the material contraction /3 and the 
rotation 0q. In Fig. [Hwe show best fits of n as function 
of #0 an d for four values of f3. 

In summary, we pointed out the role of residual stresses 
in determining the crack patterns in drying thin sub- 
strates. For spiral patterns we related the properties of 
the spiral to the degree of residual shear stress left in 
the layer. For hierarchical patterns we determined the 
position where new cracks initiate as a function of the 
mother cell, and offered a relation of final mean size of 
cells to the critical value of the yield stress. 
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